Interaction networks with significance testing

suppressMessages(suppressWarnings({
  library(magrittr, quietly = TRUE)
  library(tibble, quietly = TRUE)
  library(dplyr, quietly = TRUE)
  #library(Matrix, quietly = TRUE)
  #library(qgraph, quietly = TRUE)
  #library(psych, quietly = TRUE)
  #library(corpcor, quietly = TRUE)
  #library(corrplot, quietly = TRUE)
  library(knitr, quietly = TRUE)
}))
source('plotting_functions.R')
source('../get_cor_genes.R')
opts_chunk$set(
  dev = c("png", "pdf"), 
  dpi = 150
  #,
  #fig.width = 8/2.54, fig.height = 8/2.54
)

1 Get data

require(SummarizedExperiment)
require(magrittr)
#
#
# Join Rdata files to singe SummarizedExperiment class
#setwd('/Users/nataliabulgakova/Documents/R/RNAseq')
dir('../data', full.names = TRUE, pattern = 'Rdata') %>% sapply(function(x) get(load(x)) ) -> SE_list
SE <- do.call(cbind, SE_list)
colnames(SE) <- gsub('(SRR[0-9]+)\\.bam', '\\1', colnames(SE))

2 Hi correlations with given gene

library(org.Hs.eg.db)
require(DESeq2)
sig <- getHiCor(SE, SYMBOL='CDH1', rpkm='TRUE', cor_cutoff=0.7, anticor_cutoff=-0.6)
sig

3 Correlations

library(psych, quietly = TRUE)
M <- t(fpkm(DESeqDataSet(SE, ~1), robust = TRUE)[sig$entrez,])
colnames(M) <- sig$SYMBOL
R <- cor(M)
P  <- corr.p(R, ncol(M), adjust="fdr")$p
if(!all( R[lower.tri(R)] == t(R)[lower.tri(R)] )) stop('Matrix R is not symetric!')
P[lower.tri(P)] <- t(P)[lower.tri(P)]

3.1 Plot

library(corrplot, quietly = TRUE)
corrplot(
  R, order = "AOE", type = 'lower', diag=T, 
  method = 'number', tl.pos = 'lt', p.mat=P, sig.level=0.05
)
corrplot(
  R, order = "AOE", type = 'upper',diag=F, 
  add = TRUE, method = 'circle', tl.pos = 'n', p.mat=P, insig = "p-value", sig.level=0.05
)

3.2 Speraman correlation coeficiemnt

knitr::kable(R, "html", table.attr = 'class="table table-condensed corrtable"', digits =2, padding = 1)

	CDH1	ESRP2	MARVELD3	TSPAN13	EPCAM	GRHL2	SNRPD3	TJP3	NECTIN4	ARHGEF16	LAD1	SPINT1	TMEM141	SPRR2A	SH3YL1	PEX13	TMPRSS4	TRAF4	DHRS13	MISP	LIX1L	CCDC88A	CD99L2
CDH1	1.00	0.88	0.87	0.77	0.75	0.74	0.73	0.73	0.72	0.72	0.71	0.71	0.71	0.71	0.71	0.71	0.70	0.70	0.70	0.70	-0.65	-0.61	-0.60
ESRP2	0.88	1.00	0.83	0.68	0.73	0.84	0.58	0.75	0.73	0.81	0.62	0.79	0.65	0.50	0.68	0.57	0.52	0.72	0.56	0.69	-0.70	-0.61	-0.64
MARVELD3	0.87	0.83	1.00	0.77	0.71	0.75	0.69	0.62	0.80	0.75	0.85	0.73	0.70	0.70	0.78	0.59	0.71	0.67	0.67	0.55	-0.65	-0.63	-0.64
TSPAN13	0.77	0.68	0.77	1.00	0.59	0.46	0.66	0.56	0.76	0.64	0.70	0.48	0.72	0.80	0.69	0.72	0.69	0.60	0.70	0.50	-0.34	-0.29	-0.40
EPCAM	0.75	0.73	0.71	0.59	1.00	0.56	0.45	0.67	0.59	0.79	0.42	0.60	0.64	0.64	0.62	0.69	0.65	0.49	0.59	0.68	-0.58	-0.56	-0.60
GRHL2	0.74	0.84	0.75	0.46	0.56	1.00	0.45	0.54	0.66	0.60	0.62	0.69	0.37	0.37	0.68	0.35	0.42	0.51	0.37	0.43	-0.71	-0.65	-0.60
SNRPD3	0.73	0.58	0.69	0.66	0.45	0.45	1.00	0.40	0.53	0.43	0.64	0.48	0.65	0.54	0.42	0.57	0.52	0.47	0.62	0.41	-0.33	-0.44	-0.48
TJP3	0.73	0.75	0.62	0.56	0.67	0.54	0.40	1.00	0.69	0.74	0.44	0.63	0.58	0.57	0.59	0.54	0.62	0.70	0.56	0.75	-0.47	-0.50	-0.47
NECTIN4	0.72	0.73	0.80	0.76	0.59	0.66	0.53	0.69	1.00	0.68	0.79	0.60	0.73	0.79	0.78	0.64	0.78	0.61	0.67	0.46	-0.42	-0.48	-0.45
ARHGEF16	0.72	0.81	0.75	0.64	0.79	0.60	0.43	0.74	0.68	1.00	0.56	0.77	0.73	0.59	0.61	0.50	0.65	0.75	0.63	0.74	-0.61	-0.50	-0.55
LAD1	0.71	0.62	0.85	0.70	0.42	0.62	0.64	0.44	0.79	0.56	1.00	0.52	0.67	0.74	0.77	0.52	0.71	0.57	0.68	0.37	-0.48	-0.47	-0.47
SPINT1	0.71	0.79	0.73	0.48	0.60	0.69	0.48	0.63	0.60	0.77	0.52	1.00	0.60	0.32	0.47	0.35	0.49	0.65	0.54	0.64	-0.63	-0.64	-0.62
TMEM141	0.71	0.65	0.70	0.72	0.64	0.37	0.65	0.58	0.73	0.73	0.67	0.60	1.00	0.76	0.58	0.69	0.72	0.62	0.79	0.64	-0.27	-0.43	-0.41
SPRR2A	0.71	0.50	0.70	0.80	0.64	0.37	0.54	0.57	0.79	0.59	0.74	0.32	0.76	1.00	0.78	0.82	0.91	0.50	0.75	0.46	-0.28	-0.33	-0.32
SH3YL1	0.71	0.68	0.78	0.69	0.62	0.68	0.42	0.59	0.78	0.61	0.77	0.47	0.58	0.78	1.00	0.65	0.66	0.38	0.47	0.39	-0.56	-0.47	-0.52
PEX13	0.71	0.57	0.59	0.72	0.69	0.35	0.57	0.54	0.64	0.50	0.52	0.35	0.69	0.82	0.65	1.00	0.77	0.40	0.62	0.49	-0.31	-0.38	-0.42
TMPRSS4	0.70	0.52	0.71	0.69	0.65	0.42	0.52	0.62	0.78	0.65	0.71	0.49	0.72	0.91	0.66	0.77	1.00	0.61	0.81	0.52	-0.33	-0.43	-0.40
TRAF4	0.70	0.72	0.67	0.60	0.49	0.51	0.47	0.70	0.61	0.75	0.57	0.65	0.62	0.50	0.38	0.40	0.61	1.00	0.71	0.69	-0.44	-0.41	-0.50
DHRS13	0.70	0.56	0.67	0.70	0.59	0.37	0.62	0.56	0.67	0.63	0.68	0.54	0.79	0.75	0.47	0.62	0.81	0.71	1.00	0.52	-0.26	-0.33	-0.36
MISP	0.70	0.69	0.55	0.50	0.68	0.43	0.41	0.75	0.46	0.74	0.37	0.64	0.64	0.46	0.39	0.49	0.52	0.69	0.52	1.00	-0.48	-0.58	-0.50
LIX1L	-0.65	-0.70	-0.65	-0.34	-0.58	-0.71	-0.33	-0.47	-0.42	-0.61	-0.48	-0.63	-0.27	-0.28	-0.56	-0.31	-0.33	-0.44	-0.26	-0.48	1.00	0.60	0.64
CCDC88A	-0.61	-0.61	-0.63	-0.29	-0.56	-0.65	-0.44	-0.50	-0.48	-0.50	-0.47	-0.64	-0.43	-0.33	-0.47	-0.38	-0.43	-0.41	-0.33	-0.58	0.60	1.00	0.59
CD99L2	-0.60	-0.64	-0.64	-0.40	-0.60	-0.60	-0.48	-0.47	-0.45	-0.55	-0.47	-0.62	-0.41	-0.32	-0.52	-0.42	-0.40	-0.50	-0.36	-0.50	0.64	0.59	1.00

3.3 Cor test p-val

knitr::kable(P, "html", table.attr = 'class="table table-condensed corrtable"', digits =2, padding = 1)

	CDH1	ESRP2	MARVELD3	TSPAN13	EPCAM	GRHL2	SNRPD3	TJP3	NECTIN4	ARHGEF16	LAD1	SPINT1	TMEM141	SPRR2A	SH3YL1	PEX13	TMPRSS4	TRAF4	DHRS13	MISP	LIX1L	CCDC88A	CD99L2
CDH1	0	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
ESRP2	0	0.00	0.00	0.00	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.02	0.00	0.01	0.02	0.00	0.01	0.00	0.00	0.00	0.00
MARVELD3	0	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.01	0.00	0.00	0.00	0.01	0.00	0.00	0.00
TSPAN13	0	0.00	0.00	0.00	0.01	0.03	0.00	0.01	0.00	0.00	0.00	0.03	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.02	0.12	0.18	0.07
EPCAM	0	0.00	0.00	0.01	0.00	0.01	0.04	0.00	0.00	0.00	0.05	0.00	0.00	0.00	0.00	0.00	0.00	0.02	0.01	0.00	0.01	0.01	0.00
GRHL2	0	0.00	0.00	0.03	0.01	0.00	0.04	0.01	0.00	0.00	0.00	0.00	0.09	0.09	0.00	0.11	0.05	0.02	0.09	0.04	0.00	0.00	0.00
SNRPD3	0	0.01	0.00	0.00	0.04	0.04	0.00	0.07	0.01	0.05	0.00	0.03	0.00	0.01	0.05	0.01	0.02	0.03	0.00	0.06	0.13	0.04	0.03
TJP3	0	0.00	0.00	0.01	0.00	0.01	0.07	0.00	0.00	0.00	0.04	0.00	0.01	0.01	0.01	0.01	0.00	0.00	0.01	0.00	0.03	0.02	0.03
NECTIN4	0	0.00	0.00	0.00	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.03	0.06	0.03	0.04
ARHGEF16	0	0.00	0.00	0.00	0.00	0.00	0.05	0.00	0.00	0.00	0.01	0.00	0.00	0.01	0.00	0.02	0.00	0.00	0.00	0.00	0.00	0.02	0.01
LAD1	0	0.00	0.00	0.00	0.05	0.00	0.00	0.04	0.00	0.01	0.00	0.02	0.00	0.00	0.00	0.02	0.00	0.01	0.00	0.09	0.03	0.03	0.03
SPINT1	0	0.00	0.00	0.03	0.00	0.00	0.03	0.00	0.00	0.00	0.02	0.00	0.00	0.13	0.03	0.11	0.02	0.00	0.01	0.00	0.00	0.00	0.00
TMEM141	0	0.00	0.00	0.00	0.00	0.09	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.21	0.05	0.06
SPRR2A	0	0.02	0.00	0.00	0.00	0.09	0.01	0.01	0.00	0.01	0.00	0.13	0.00	0.00	0.00	0.00	0.00	0.02	0.00	0.03	0.20	0.12	0.14
SH3YL1	0	0.00	0.00	0.00	0.00	0.00	0.05	0.01	0.00	0.00	0.00	0.03	0.01	0.00	0.00	0.00	0.00	0.08	0.03	0.08	0.01	0.03	0.02
PEX13	0	0.01	0.01	0.00	0.00	0.11	0.01	0.01	0.00	0.02	0.02	0.11	0.00	0.00	0.00	0.00	0.00	0.07	0.00	0.02	0.15	0.08	0.05
TMPRSS4	0	0.02	0.00	0.00	0.00	0.05	0.02	0.00	0.00	0.00	0.00	0.02	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.02	0.12	0.05	0.07
TRAF4	0	0.00	0.00	0.00	0.02	0.02	0.03	0.00	0.00	0.00	0.01	0.00	0.00	0.02	0.08	0.07	0.00	0.00	0.00	0.00	0.04	0.06	0.02
DHRS13	0	0.01	0.00	0.00	0.01	0.09	0.00	0.01	0.00	0.00	0.00	0.01	0.00	0.00	0.03	0.00	0.00	0.00	0.00	0.02	0.22	0.13	0.10
MISP	0	0.00	0.01	0.02	0.00	0.04	0.06	0.00	0.03	0.00	0.09	0.00	0.00	0.03	0.08	0.02	0.02	0.00	0.02	0.00	0.03	0.01	0.02
LIX1L	0	0.00	0.00	0.12	0.01	0.00	0.13	0.03	0.06	0.00	0.03	0.00	0.21	0.20	0.01	0.15	0.12	0.04	0.22	0.03	0.00	0.00	0.00
CCDC88A	0	0.00	0.00	0.18	0.01	0.00	0.04	0.02	0.03	0.02	0.03	0.00	0.05	0.12	0.03	0.08	0.05	0.06	0.13	0.01	0.00	0.00	0.01
CD99L2	0	0.00	0.00	0.07	0.00	0.00	0.03	0.03	0.04	0.01	0.03	0.00	0.06	0.14	0.02	0.05	0.07	0.02	0.10	0.02	0.00	0.01	0.00

4 Speraman correlation network - spring layout

4.1 Significant correlations

library(qgraph, quietly = TRUE)
  r <- R
  r[P > 0.05] <- NA 
  diag(r) <- NA
  z <- qgraph(r, layout = "spring", title = 'Speraman correlation network - significant correlations', DoNotPlot=TRUE, directed=FALSE, bidirectional=TRUE)
  z$graphAttributes$Nodes$labels <- names(z$graphAttributes$Nodes$labels)
  plot(z)

4.2 All correlations

  z <- qgraph(R, layout = "spring", title = 'Speraman correlation network - all correlations', DoNotPlot=TRUE, directed=FALSE, bidirectional=TRUE)
  z$graphAttributes$Nodes$labels <- names(z$graphAttributes$Nodes$labels)
  plot(z)

5 Speraman correlation network - circle layout

5.1 Significant correlations

  z <- qgraph(r, layout = "circle", title = 'Speraman correlation network - significant correlations', DoNotPlot=TRUE, directed=FALSE, bidirectional=TRUE)
  z$graphAttributes$Nodes$labels <- names(z$graphAttributes$Nodes$labels)
  plot(z)

5.2 All correlations

  z <- qgraph(R, layout = "circle", title = 'Speraman correlation network - all correlations', DoNotPlot=TRUE, directed=FALSE, bidirectional=TRUE)
  z$graphAttributes$Nodes$labels <- names(z$graphAttributes$Nodes$labels)
  plot(z)

6 Speraman correlation network - with LASSO

if(matrixcalc::is.positive.definite(R)) {
  r <- R
  p <- P
} else {
    #Needed only if cor(M) is no PD
  colnames(M) <- sig$SYMBOL
  MM <- collapsemaxcor(M, method = 'pearson')
  while(!attr(MM, 'PD')) {
    MM <- collapsemaxcor(MM, method = 'pearson')
  }
  
  r <- cor(M)
  p <- corr.p(r, 20000, adjust = 'none')$p
}
  

6.1 Significant correlations

  r[p > 0.05] <- 0
 
  
  if(matrixcalc::is.positive.definite(R)) {
    r <- as.matrix(nearPD(r)$mat)
  }
 EBICgraph <- EBICglasso(r, 20000, 0.5)
  
  z <- qgraph(EBICgraph, layout = "spring", title = 'Speraman correlation network - significant correlations', DoNotPlot=TRUE, directed=FALSE, bidirectional=TRUE)
  z$graphAttributes$Nodes$labels <- names(z$graphAttributes$Nodes$labels)
  plot(z)

6.2 All correlations

  r <- cor(M)
  EBICgraph <- EBICglasso(r, 20000, 0.5)

Network with lowest lambda selected as best network. Try setting 'lambda.min.ratio' lower.

  z <- qgraph(EBICgraph, layout = "spring", title = 'Speraman correlation network - all correlations', DoNotPlot=TRUE, directed=FALSE, bidirectional=TRUE)
  z$graphAttributes$Nodes$labels <- names(z$graphAttributes$Nodes$labels)
  plot(z)